//Vector_4
template <typename scalar>
void Vector_4<scalar> :: Invert()
{
  x *= -1;
  y *= -1;
  z *= -1;
  w *= -1;
}

template <typename scalar>
void Vector_4<scalar> ::Normalize()
{
	scalar l = sqrt (com::sqr(x) + com::sqr(y) + com::sqr(z) + com::sqr(w));
	x /= l;
	y /= l;
	z /= l;
	w /= l;
}
//internal operators
template <typename scalar>
scalar Vector_4<scalar> :: operator [](size_t idx) const 
{
	switch (idx)
	{
	case 0: 
    {
      return x; 
      break;
    }
  case 1: 
    {
      return y; 
      break;
    }
  case 2: 
    {
      return z; 
      break;
    }
  case 3:
    {
      return w;
      break;
    }
  default:
    return 0;
	};
}

template <typename scalar>
Vector_4<scalar>& Vector_4<scalar> :: operator += (const Vector_4<scalar>& p)
{
  x += p.x;
  y += p.y;
  z += p.z;
  w += p.w;
  return *this;
}

template <typename scalar>
Vector_4<scalar>& Vector_4<scalar> :: operator -= (const Vector_4<scalar>& p)
{
  x -= p.x;
  y -= p.y;
  z -= p.z;
  w -= p.w;
  return *this;
}

template <typename scalar>
Vector_4<scalar>& Vector_4<scalar> :: operator /= (const scalar& k)
{
  x /= k;
  y /= k;
  z /= k;
  w /= k;
  return *this;
}

template <typename scalar>
Vector_4<scalar>& Vector_4<scalar> :: operator %= (const Vector_4<scalar>& p)
{
	return Vector_4<scalar>(x com::mod p.x, y com::mod p.y, z com::mod p.z, w com::mod p.w);
}

template <typename scalar>
Vector_4<scalar>& Vector_4<scalar> :: operator &= (const Vector_4<scalar>& p)
{
	return Vector_4<scalar>(x * p.x, y * p.y, z * p.z, w * p.w);
}

template <typename scalar>
mBool Vector_4<scalar> :: operator != (const Vector_4<scalar>& b)
{
	if (x == b.x && y == b.y && z == b.z && w == b.w)
	{return false;}
	else
	{return true;};
}

template <typename scalar>
mBool Vector_4<scalar> :: operator == (const Vector_4<scalar>& b)
{
	if (x == b.x && y == b.y && z == b.z && w == b.w)
	{return true;}
	else
	{return false;};
}

//external operators

template <typename scalar>
const Vector_4<scalar> operator + (const Vector_4<scalar>& a, const Vector_4<scalar>& b) 
{
	return Vector_4<scalar>(a.x + b.x, a.y + b.y, a.z + b.z, a.w + b.w);
}

template <typename scalar>
const Vector_4<scalar> operator - (const Vector_4<scalar>& a, const Vector_4<scalar>& b)
{
	return Vector_4<scalar>(a.x - b.x, a.y - b.y, a.z - b.z, a.w - b.w);
}

template <typename scalar>
const Vector_4<scalar> operator / (const Vector_4<scalar>& a, const scalar& k)
{
	return Vector_4<scalar>(a.x / k, a.y / k, a.z / k, a.w / k);
}

template <typename scalar>
const Vector_4<scalar> operator ^ (const Vector_4<scalar>& a, const Vector_4<scalar>& b)
{
	return Vector_4<scalar>(a.y * b.z - b.y * a.z, a.z * b.x - b.z * a.x, a.x * b.y - b.x * a.y, scalar(0));
}

template <typename scalar> 
const scalar operator * (const Vector_4<scalar>& a, const Vector_4<scalar>& b)
{
	return a.x * b.x + a.y * b.y + a.z * b.z + a.w * b.w;
}

template <typename scalar>
const Vector_4<scalar> operator * (const Vector_4<scalar>& a, const scalar& k)
{
	return Vector_4<scalar>(a.x * k, a.y * k, a.z * k, a.w * k);
}

template <typename scalar>
const Vector_4<scalar> operator * (const scalar& k, const Vector_4<scalar>& a)
{
	return Vector_4<scalar>(a.x * k, a.y * k, a.z * k, a.w * k);
}

template <typename scalar>
const Vector_4<scalar> operator % (const Vector_4<scalar>& a, const Vector_4<scalar>& b)
{
	return Vector_4<scalar>(a.x com::mod b.x, a.y com::mod b.y, a.z com::mod b.z, a.w com::mod b.w);
}

template <typename scalar>
const Vector_4<scalar> operator & (const Vector_4<scalar>& a, const Vector_4<scalar>& b)
{
	return Vector_4<scalar>(a.x * b.x, a.y * b.y, a.z * b.z, a.w * b.w);
}

//Vector_3
//internal operators

template <typename scalar> 
Vector_3<scalar>& Vector_3<scalar> :: operator += (const Vector_3<scalar>& p)
{
  x += p.x;
  y += p.y;
  z += p.z;
  return *this;
}

template <typename scalar>
Vector_3<scalar>& Vector_3<scalar> :: operator -= (const Vector_3<scalar>& p)
{
	x -= p.x;
	y -= p.y;
	z -= p.z;
  return *this;
}


template <typename scalar>
Vector_3<scalar>& Vector_3<scalar> :: operator /= (const scalar& k)
{
	x /= k;
	y /= k;
	z /= k;
	return *this;
}

template <typename scalar>
Vector_3<scalar>& Vector_3<scalar> :: operator ^= (const Vector_3<scalar>& p)
{
  return Vector_3<scalar>(y * p.z - p.y * z, z * p.x - p.z * x, x * p.y - p.x * y);
}


template <typename scalar>
Vector_3<scalar>& Vector_3<scalar> :: operator %= (const Vector_3<scalar>& p)
{
	return Vector_3<scalar>(x com::mod p.x, y com::mod p.y, z com::mod p.z);
}

template <typename scalar>
Vector_3<scalar>& Vector_3<scalar> :: operator &= (const Vector_3<scalar>& p)
{
	return Vector_3<scalar>(x * p.x, y * p.y, z * p.z);
}

template <typename scalar>
mBool Vector_3<scalar> :: operator != (const Vector_3<scalar>& b)
{
	if (x == b.x && y == b.y && z == b.z)
	{return false;}
	else
	{return true;};
}

template <typename scalar>
mBool Vector_3<scalar> :: operator == (const Vector_3<scalar>& b)
{
	if (x == b.x && y == b.y && z == b.z)
	{return true;}
	else
	{return false;};
}

template <typename scalar>
Vector_3<scalar>& Vector_3<scalar> :: operator *= (const scalar& k)
{
	x *= k;
	y *= k;
	z *= k;
	return *this;
}

//external operators

template <typename scalar>
const Vector_3<scalar> operator + (const Vector_3<scalar>& a, const Vector_3<scalar>& b) 
{
    return Vector_3<scalar>(a.x + b.x, a.y + b.y, a.z + b.z);
}

template <typename scalar>
const Vector_3<scalar> operator - (const Vector_3<scalar>& a, const Vector_3<scalar>& b)
{
	return Vector_3<scalar>(a.x - b.x, a.y - b.y, a.z - b.z);
}

template <typename scalar>
const Vector_3<scalar> operator / (const Vector_3<scalar>& a, const scalar& k)
{
	return Vector_3<scalar>(a.x / k, a.y / k, a.z / k);
}

template <typename scalar>
const Vector_3<scalar> operator ^ (const Vector_3<scalar>& a, const Vector_3<scalar>& b)
{
	return Vector_3<scalar>(a.y * b.z - b.y * a.z, a.z * b.x - b.z * a.x, a.x * b.y - b.x * a.y);
}

template <typename scalar> 
const scalar operator * (const Vector_3<scalar>& a, const Vector_3<scalar>& b)
{
    return a.x * b.x + a.y * b.y + a.z * b.z;
}

template <typename scalar>
const Vector_3<scalar> operator * (const Vector_3<scalar>& a, const scalar& k)
{
	return Vector_3<scalar>(a.x * k, a.y * k, a.z * k);
}

template <typename scalar>
const Vector_3<scalar> operator * (const scalar& k, const Vector_3<scalar>& a)
{
	return Vector_3<scalar>(a.x * k, a.y * k, a.z * k);
}

template <typename scalar>
const Vector_3<scalar> operator % (const Vector_3<scalar>& a, const Vector_3<scalar>& b)
{
	return Vector_3<scalar>(a.x com::mod b.x, a.y com::mod b.y, a.z com::mod b.z);
}

template <typename scalar>
const Vector_3<scalar> operator & (const Vector_3<scalar>& a, const Vector_3<scalar>& b)
{
	return Vector_3<scalar>(a.x * b.x, a.y * b.y, a.z * b.z);
}

template <typename scalar>
scalar& Vector_3<scalar> :: operator [](size_t idx)
{
	switch (idx)
	{
	case 0: 
    {
      return x; 
      break;
    };
	case 1: 
    {
    return y; 
    break;
    };
	case 2: 
    {
      return z; 
      break;
    };
  default:
    break;
	};
  return z;
}

template <typename scalar>
scalar Vector_3<scalar> :: operator [](size_t idx) const
{
	switch (idx)
	{
	case 0: 
    {
      return x; 
      break;
    };
	case 1: 
    {
    return y; 
    break;
    };
	case 2: 
    {
      return z; 
      break;
    };
  default:
    return 0;
    break;
	};
}

//functions
template <typename scalar>
const Vector_3<scalar> Interpolate (const Vector_3<scalar>& a, const Vector_3<scalar>& b, const mFloat& k)
{
	return Vector_3<scalar>(a.x + (b.x - a.x) / k, a.y + (b.y - a.y) / k, a.z + (b.z - a.z) / k);
}

template <typename scalar>
void Vector_3<scalar> :: Invert (void)
{
  x *= -1;
  y *= -1;
  z *= -1;
}

template <typename scalar>
void Vector_3<scalar> :: Normalize (void)
{
	scalar l = sqrt (com::sqr(x) + com::sqr(y) + com::sqr(z));
  x /= l;
	y /= l;
	z /= l;
}

template <typename scalar>
void Vector_3<scalar> :: SetLength (scalar& length)
{
	scalar old = sqrt (com::sqr(x) + com::sqr(y) + com::sqr(z));
	x /= old;
	y /= old;
	z /= old;
	x *= length;
  y *= length;
	z *= length;
}

template <typename scalar>
void Vector_3<scalar> :: RotateAround (const Vector_3<scalar>& target, const mFloat& angle)
{
	
  Vector_3<scalar> v = *this;
  scalar x_, y_, z_, temp1;
  scalar rad_angle = com::grad2rad(angle);
  scalar c = cos(rad_angle);
  scalar s = sin(rad_angle);

  temp1 = sqrt(1 - com::sqr(target.z));
  if (temp1 != 0)
  {
    x_ = v.x * target.y / temp1 - v.y * target.x / temp1;
    y_ = v.x * target.x / temp1 + v.y * target.y / temp1;
    v.x = x_;
    v.y = y_;
    y_ = v.y * target.z - v.z * temp1;
    z_ = v.y * temp1 + v.z * target.z;
    v.y = y_;
    v.z = z_;
    x_ = c * v.x - s * v.y;
    y_ = c * v.y + s * v.x;
    v.x = x_;
    v.y = y_;
    y_ = v.y * target.z + v.z * temp1;
    z_ = -v.y * temp1 + v.z * target.z;
    v.y = y_;
    v.z = z_;
    x_ = v.x * target.y / temp1 + v.y * target.x / temp1;
    y_ = -v.x * target.x / temp1 + v.y * target.y / temp1;
    v.x = x_;
    v.y = y_;
    v.Normalize();
  }
  else
  {
    x_ = c * v.x - s * v.y;
    y_ = c * v.y + s * v.x;
    v.x = x_;
    v.y = y_;
  };
  *this = v;
}


//Vector_2<scalar>


//internal operators

template <typename scalar> 
Vector_2<scalar>& Vector_2<scalar> :: operator += (const Vector_2<scalar>& p)
{
    x += p.x;
    y += p.y;
    return *this;
}

template <typename scalar>
Vector_2<scalar>& Vector_2<scalar> :: operator -= (const Vector_2<scalar>& p)
{
	x -= p.x;
	y -= p.y;
  return *this;
}


template <typename scalar>
Vector_2<scalar>& Vector_2<scalar> :: operator /= (const scalar& k)
{
	x /= k;
	y /= k;
	return *this;
}



template <typename scalar>
Vector_2<scalar>& Vector_2<scalar> :: operator %= (const Vector_2<scalar>& p)
{
	return Vector_2<scalar>(x com::mod p.x, y com::mod p.y);
}

template <typename scalar>
Vector_2<scalar>& Vector_2<scalar> :: operator &= (const Vector_2<scalar>& p)
{
	return Vector_2<scalar>(x * p.x, y * p.y);
}

template <typename scalar>
mBool Vector_2<scalar> :: operator != (const Vector_2<scalar>& b)
{
	if (x == b.x && y == b.y)
	{return false;}
	else
	{return true;};
}

template <typename scalar>
mBool Vector_2<scalar> :: operator == (const Vector_2<scalar>& b)
{
	if (x == b.x && y == b.y)
	{return true;}
	else
	{return false;};
}

template <typename scalar>
Vector_2<scalar>& Vector_2<scalar> :: operator *= (const scalar& k)
{
	x *= k;
	y *= k;
	return *this;
}

//external operators

template <typename scalar>
const Vector_2<scalar> operator + (const Vector_2<scalar>& a, const Vector_2<scalar>& b) 
{
    return Vector_2<scalar>(a.x + b.x, a.y + b.y);
}

template <typename scalar>
const Vector_2<scalar> operator - (const Vector_2<scalar>& a, const Vector_2<scalar>& b)
{
	return Vector_2<scalar>(a.x - b.x, a.y - b.y);
}

template <typename scalar>
const Vector_2<scalar> operator / (const Vector_2<scalar>& a, const scalar& k)
{
	return Vector_2<scalar>(a.x / k, a.y / k);
}


template <typename scalar> 
const scalar operator * (const Vector_2<scalar>& a, const Vector_2<scalar>& b)
{
    return a.x * b.x + a.y * b.y;
}

template <typename scalar>
const Vector_2<scalar> operator * (const Vector_2<scalar>& a, const scalar& k)
{
	return Vector_2<scalar>(a.x * k, a.y * k);
}

template <typename scalar>
const Vector_2<scalar> operator * (const scalar& k, const Vector_2<scalar>& a)
{
	return Vector_2<scalar>(a.x * k, a.y * k);
}

template <typename scalar>
const Vector_2<scalar> operator % (const Vector_2<scalar>& a, const Vector_2<scalar>& b)
{
	return Vector_2<scalar>(a.x com::mod b.x, a.y com::mod b.y);
}

template <typename scalar>
const Vector_2<scalar> operator & (const Vector_2<scalar>& a, const Vector_2<scalar>& b)
{
	return Vector_3<scalar>(a.x * b.x, a.y * b.y);
}

template <typename scalar>
scalar Vector_2<scalar> :: operator [](size_t idx) const 
{
	switch (idx)
	{
	case 0: {return x; break;};
	case 1: {return y; break;};
	};
}

//functions
template <typename scalar>
const Vector_2<scalar> Interpolate (const Vector_2<scalar>& a, const Vector_2<scalar>& b, const mFloat& k)
{
	return Vector_2<scalar>(a.x + (b.x - a.x) / k, a.y + (b.y - a.y) / k);
}

template <typename scalar>
void Vector_2<scalar> :: Invert (void)
{
    x *= -1;
	y *= -1;
}

template <typename scalar>
void Vector_2<scalar> :: Normalize (void)
{
	scalar l = sqrt (sqr(x) + sqr(y));
    x /= l;
	y /= l;
}

template <typename scalar>
void Vector_2<scalar> :: SetLength (scalar& length)
{
	scalar old = sqrt (sqr(x) + sqr(y));
	x /= old;
	y /= old;
	
	x *= length;
    y *= length;
	
}

template <typename scalar>
void Vector_2<scalar> :: Rotate(const mFloat& angle)
{
	
  mFloat X_, Y_, Z_;
  mFloat rad_angle = com::grad2rad(angle);
  mFloat s = sin(rad_angle);
  mFloat c = cos(rad_angle);
    //=================== rotate XY ==============================
    X_ = X * c - Y * s;
    Y_ = Y * c + X * s;
    x = X_;
    y = Y_;
    Normalize();
}